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Solitary waves bifurcated from edges of Bloch bands in two-dimensional periodic media are deter- 
mined both analytically and numerically in the context of a two-dimensional nonlinear Schrodinger 
equation with a periodic potential. Using multi-scale perturbation methods, envelope equations of 
solitary waves near Bloch bands are analytically derived. These envelope equations reveal that soli- 
tary waves can bifurcate from edges of Bloch bands under either focusing or defocusing nonlinearity, 
depending on the signs of second-order dispersion coefficients at the edge points. Interestingly, at 
] edge points with two linearly independent Bloch modes, the envelope equations lead to a host of 

, solitary wave structures including reduced-symmetry solitons, dipole-array solitons, vortex-cell soli- 

■ tons, and so on — many of which have never been reported before. It is also shown analytically that 

£\j | the centers of envelope solutions can only be positioned at four possible locations at or between po- 

tential peaks. Numerically, families of these solitary waves are directly computed both near and far 
away from band edges. Near the band edges, the numerical solutions spread over many lattice sites, 
and they fully agree with the analytical solutions obtained from envelope equations. Far away from 
the band edges, solitary waves are strongly localized with intensity and phase profiles characteristic 
of individual families. 
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Qh I I. INTRODUCTION 

Nonlinear wave propagation in periodic media is attracting a lot of attention these days. This was stimulated 
in part by rapid advances in optics, Bose-Einstein condensates, and related fields. In optics, various periodic and 
quasi-periodic structures (such as photonic crystals, photonic crystal fibers, periodic waveguide arrays and photonic 
lattices) have been constructed by ingenious experimental techniques, with applications to light routing, switching 
' and optical information processing [J S S El IS B 0> 111 • These periodic media create a wide range of new phenomena 
for light propagation, even in the linear regime. For instance, the diffraction of light in a periodic medium exhibits 
distinctively different patterns from homogeneous diffraction If the periodic medium has a local defect, this 
defect can guide light by a totally new physical mechanism called repeated Bragg reflections [l|, 0, [t| [l(| [HI US, EH • 
When the nonlinear effects become significant, such as with high-power beams or strongly nonlinear materials, the 
physical phenomena would be even richer and more complex, and their understanding is far from complete yet. In 
i Bose-Einstein condensates, one direction of recent research is to load the condensates into periodic optical lattices 
[Til [l5l Hil . This problem and the above nonlinear optics problems are closely related, and are often analyzed together. 

Solitary waves play an important role in nonlinear wave systems. These waves are nonlinear localized structures 
which propagate without change of shape. In physical communities, they are often just called solitons, which we do 
occasionally in this paper as well. In one-dimensional (ID) periodic media, solitary waves (called lattice solitons) 
have been predicted and observed in optical experiments [jj, H, E3, EM 0, l20| . But in two and higher dimensions, 
periodic media can support a much wider array of solitary wave structures, many of which have no counterparts in 
ID systems. So far, several types of 2D lattice solitons in the semi-infinite bandgap (such as fundamental, dipole 
and vortex solitons) as well as the first bandgap (such as fundamental, vortex and reduced-symmetry gap solitons) 
have been reported % El [H, [H, H, EEM H [H, HI HI . Solitons in Bessel-ring lattices and 2D quasi- 
periodic lattices have been reported as well [32l. l33l. l34j ■ All these works were either numerical or experimental, and an 
analytical understanding on these solitons is still lacking. Some of these solitons bifurcate from edges of Bloch bands. 
For instance, the gap vortex solitons reported in [27| bifurcate from the two X-symmetry points of the second Bloch 
band, and the reduced-symmetry solitons reported in [301 ] bifurcate from a single X-symmetry point of the second 
Bloch band (both under focusing nonlinearity). This raises the following important questions: are there other types 
of solitary waves bifurcated from Bloch bands in 2D periodic media? how can such solitary waves be analytically 
predicted and classified? 

In this paper, we determine all possible solitary-wave structures bifurcated from edges of Bloch bands in 2D periodic 
media both analytically and numerically, using the two-dimensional nonlinear Schrodinger equation with a periodic 
potential as the mathematical model. By multi-scale perturbation methods, we derive the envelope equations of 
these solitary waves near band edges. We find that these envelope equations admit solutions which lead to not only 
solitons reported before (see [27l. |30| for instance), but also new solitary- wave structures such as dipole-array solitons 
bifurcated from the second Bloch band under focusing nonlinearity, and vortex-cell solitons bifurcated from the second 
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Bloch band under defocusing nonlinearity. We also show analytically that the centers of envelope solutions can only 
be positioned at four possible locations at or between the potential peaks. We further determine the whole families of 
these solitons both near and far away from band edges directly using numerical methods f35j. Near the band edges, 
the numerical solutions spread over many lattice sites, and they fully agree with the analytical solutions obtained 
from envelope equations. Far away from the band edges, solitary waves are strongly localized, and their intensity and 
phase profiles carry signatures of individual soliton families. These studies provide a rather complete understanding 
of solitary waves bifurcated from Bloch bands in 2D periodic media. 

II. THE MATHEMATICAL MODEL 

The mathematical model we use for the study of solitary waves in 2D periodic media is the 2D nonlinear Schr" odinger 
(NLS) equation with a periodic potential: 

iU t + U xx + U yy - V(x, y)U + a\U\ 2 U = 0, (2.1) 

where U(x, y, t) is a complex function, V(x, y) is the periodic potential (also called the lattice potential), and a — ±1 
is the sign of nonlinearity. This model arises in Bose-Einstein condensates trapped in a 2D optical lattice (where t 
is time) HI [H as well as light propagation in a periodic Kerr medium under paraxial approximation (where t is 
the distance of propagation). In certain optical materials (such as photorefractive crystals), the nonlinearity is of a 
different (saturable) type. But it is known that those different types of nonlinearities give qualitatively similar results 
as the Kerr nonlinearity above [U, [U [H, [3(| ■ 
In this article we take the lattice potential as 

V(x, y) = V Q (sin 2 x + sin 2 y) , (2.2) 

whose periods L along the x and y directions are both equal to 7r. This square-lattice potential can be readily engi- 
neered in Bose-Einstein condensates [HI, HH and optics (flol • This potential is separable, which makes our theoretical 
analysis a little easier. Similar analysis can be repeated for other types of periodic potentials and nonlinearities (such 
as saturable nonlinearities in photorefractive crystals [H,[36j]) with minimal changes. Without loss of generality, when 
we carry out specific computations, we always set Vq = 6 in the potential (12. 2p . 
Solitary waves in Eq. (12. 1|) are sought in the form 

U(x,y,t)=u(x,y)e- i » t , (2.3) 

where amplitude function u(x, y) is a solution of the following equation: 

[F(x) + F(y)}u + Liu + <r\u\ 2 u = 0, (2.4) 

F(x) = Vo sin 2 x, (2.5) 

and li is a propagation constant. 

In this article, we determine solitary waves in Eq. (|2.4p which are bifurcated from the Bloch bands (i.e. continuous 
spectrum) of that equation. To do so, information about the Bloch bands of Eq. (|2.4p is essential. Such Bloch bands 
will be analyzed first below. 

III. BLOCH BANDS AND BAND GAPS 

When function u(x,y) is infinitesimal, Eq. (|2.4|) becomes a linear equation: 

u xx + Uyy - [F(x) + F(y)]u + fiu = 0. (3.1) 

Solutions of this linear equation are the Bloch modes, and the corresponding propagation constants fi form Bloch 
bands. Since the potential in (|3.ip is separable, Bloch solutions and Bloch bands of this 2D equation can be constructed 
from solutions of a ID equation. Specifically, the 2D Bloch solution u(x, y) of Eq. (|3.ip and its propagation constant 
ix can be split into the following form: 



u{x,y) = p(x;u! a )p(y;uj b ), Li = ui a +uj b , 



(3.2) 
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FIG. 1: (color online) (a) Dispersion curves of the ID equation (|3.3[) with Vb = 6; (b) Bloch bands (shaded regions) and 
bandgaps at various values of potential levels Vb in Eq. (13. 3[) . The circle points in (a) correspond to edges of Bloch bands 
marked by numbers 1-5 in (b). Bloch modes at these locations are displayed in Fig. [2] 



where p(x; to) is a solution of the following ID equation: 

Pxx ~ F(x)p + up = 0. (3.3) 
This ID equation is equivalent to the Mathieu's equation. Its solution is 

p(x;uj) = e lkx p(x;u), (3.4) 

where p(x;u>) is periodic with the same period 7r as the potential F(x), and u = uj(k) is the ID dispersion relation. 
This dispersion diagram is shown in Fig. QJa) (for Vb = 6). The bandgap structure of this ID equation (|3.3[) at 
various values of Vb is shown in Fig. [ljb). Notice that in the ID case, at any non-zero value of Vb, bandgaps appear; 
in addition, the number of bandgaps is infinite. The first five Bloch waves p(x; Uk), 1 < k < 5 at the lowest five edges 
of Bloch bands u — u>f. are displayed in Fig. [5J These Bloch waves have been normalized to have unit amplitude. 
Notice that these solutions at band edges are all real- valued. 

Using these ID dispersion results and the above connection between ID and 2D Bloch solutions, we can construct 
the dispersion surfaces and bandgap structures of the 2D problem (|3.1| . The 2D Bloch- mode solution is of the form 

u(x, y) = e lk * x+lk yyp[x- u(k x )] p[y; u(k y )}, (3.5) 

where 



H = u(k x ) + u(k y ), -l<k x ,k y <\ (3.6) 

is the 2D dispersion relation, and — 1 < k x , k y < 1 is the first Brillouin zone. This 2D dispersion relation (at Vq = 6) 
is shown in Fig. [3] (a). The 2D bandgap structure at various values of Vb is shown in Fig. [3] (b). Unlike the ID case, 
for a given Vb value, there are only a finite number of bandgaps in the 2D problem. The first bandgap appears only 
when Vb > 1.40, the second bandgap appearing when Vb > 4.13, etc. As Vb increases further, more bandgaps will be 
found. At the potential strength Vb = 6 as used in this paper, two bandgaps exist. Edges of Bloch bands at this Vb 
value are marked in Fig. [3] (b) as A, B, C, D, E' respectively. 

Locations of Bloch-band edges in the first Brillouin zone are important, as these locations reveal the symmetry 
properties of Bloch modes. To clearly mark such locations, we plotted the first Brillouin zone in Fig. [31 In the 
literature, the center of this Brillouin zone is called the T point, the four corners called the M points, and points 
(k x ,k y ) = (1,0), (0, 1) called the X and X' points respectively [27j | . These points are marked in the Brillouin zone 
of Fig. [3j Note that the four corner (M) points correspond to the same Bloch mode, but points X and X' lead to 
different (linearly independent) Bloch modes. In this first Brillouin zone, Bloch-band edges A, B in Fig. [3] (b) are 
located at T and M symmetry points respectively. At band edge C, there are two points on the dispersion surfaces 
which are located at X and X' symmetry points. The same goes to the band edge D. Band edge E is located at the 
r symmetry point. 

Now we examine 2D Bloch solutions at band edges. To illustrate, we consider the edge points A, B, C, D and E 
in Fig. [3](b), where Vb = 6. At both points A and B, there is a single Bloch solution. The Bloch solution at point 
A (located at T point of the Brillouin zone) is u(x,y) = p(x; Ui)p{y; Ui), where p{x;lu 1 ) is shown in Fig. ^MX)- This 
solution is displayed in Fig. SUA). Its propagation constant is /i = 2u\, Notice that this solution has the symmetry 
u(x,y) = u(y,x). For convenience, we denote point A as "1 + 1", and is so indicated in Fig. [3](b). Similarly, the Bloch 
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FIG. 2: One-dimensional Bloch waves of Eq. (|3.3[1 at the lowest five edges of Bloch bands marked by circles in Fig. QJa) and 
by numbers in Fig. \£[b). (1): wi = 2.063182; (2): lo 2 = 2.266735; (3): uj s = 5.165940; (4): cj 4 = 6.81429; (5) w 6 = 7.74678. 
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FIG. 3: (color online) (a) Dispersion surfaces of the 2D problem (|3.ip at Vo = 6; the low figure at its side is the first Brillouin 
zone; (b) the 2D bandgap structure for various values of Vo. Letters 'A, B, C, D, E' mark the edges of Bloch bands at Vo = 6. 
Insets like 'A:l+1' are explained in the text. 



solution at point B (located at M point of the Brillouin zone) is u(x, y) — p{x;uj2)p{y',^2), where p(x; CJ2) is shown 
in Fig. [2^2). This solution is displayed in Fig. [4f F3) , and its propagation constant is /i = 2cl>2- This solution has the 
symmetry u(x,y) = u(y,x) as well. Point B is "2 + 2" in our notations. Points C, D, E are different from A and B 
and are more interesting. At these points, there are two linearly independent Bloch solutions, u(x,y) and u(y,x). At 
point C, these two solutions are u{x,y) = p(x; u>i)p(y; u>s) and u(y,x) = p{y;uii)p{x;oj^), with the same propagation 
constant fx = u>± + W3. Here p{x\ui^,) is shown in Fig. [2^3). These solutions correspond to the X and X' points in 
the Brillouin zone. Point C is thus "1+3". Its u(x,y) solution is displayed in Fig. fJJC); the u(y,x) solution is just 
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FIG. 4: (color online) Upper left: the potential function —V(x, y); (A, B, C, D, E) are Bloch modes at the edge points of Bloch 
bands marked by the same letters in Fig. [3^b). 



a 90° rotation of u(x,y) in Fig. 0|C) and thus not shown. Point D is "2+4", where the two linearly independent 
Bloch solutions are u(x, y) = p(x; L02)p{y] ^i) and u(y, x) = p(y; u>2)p(x; u>i), the former of which is displayed in Fig. 
IDJD). These solutions also correspond to the X and X' points in the Brillouin zone. Point E is "1+5", and the two 
Bloch solutions are u(x, y) = p(x; LUi)p(y; uj^) and u(y, x) = p(y; uj\)p{x\ u>s) (the former is shown in Fig. HfE)). These 
solutions correspond to the T symmetry point in the Brillouin zone. It is interesting that at point E, the two different 
Bloch modes are both located at the same T point, while at points C and D, their two Bloch modes are located at 
different symmetry points (X and X') in the Brillouin zone. This difference between points E and {C, D} will manifest 
itself later in solitary wave bifurcations. Because of the existence of two linearly independent Bloch solutions at band 
edges C, D and E, their linear superpositions remain a solution. These superpositions can give rise to interesting 
solution patterns, some of which (such as vortex arrays) have been pointed out before [H, while many others are 
new. Solitary waves bifurcated from these linearly superimposed Bloch modes are one of the subjects we will focus 
on in the remainder of the paper. 

The above Bloch solutions exist on band edges with infinitesimal amplitudes. When amplitudes of these solutions 
increase, these Bloch solutions may localize and form solitary-wave structures. The corresponding propagation con- 
stant [i then moves from band edges into band gaps. In the next section, we analyze how solitary waves bifurcate 
from Bloch solutions of band edges using multi-scale perturbation methods. 
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IV. ENVELOPE EQUATIONS OF BLOCH MODES 



In this section, we develop an asymptotic theory to analyze small-amplitude solitary waves bifurcating from Bloch 
waves near band edges in Eq. (|2.4[) . and derive their envelope equations. We have known that at an edge of a Bloch 
band, there may be either one or two linearly independent Bloch modes. The bifurcation analysis for these two cases 
are similar, thus we will present detailed calculations only for the latter case (which is a little more complex), and 
just give the results for the former case. 



A. Derivation of envelope equations 



Let us consider a 2D band edge fio = ^o.i + ^o,2, where o>o.n (n — 1, 2) are ID band edges, uq.i ^ ^0,2, and the two 
linearly independent Bloch modes are Pi(x)p 2 (y) and Pi(y)p 2 (x) with p n (x) — p(x; uiQ >n ). Notice that 

p n (x + L) = ±p n {x) (4.1) 

since wo, n is a ID band edge. Here L — it is the period of the ID potential F(x). When the solution u(x,y) of Eq. 
(|2.4[) is infinitesimally small, this solution on the band edge is then a linear superposition of these two Bloch modes 
in the general case. When u(x,y) is small but not infinitesimal, we can expand the solution u(x,y) of (|2.4|) into a 
multi-scale perturbation series: 

u = eu + e 2 ui + <?u 2 H , (4.2) 

^i = /i + ?ye 2 , (4.3) 

where 

u = A 1 (X,Y) Pl (x)p 2 (y) + A 2 (X 7 Y)p 2 (x) Pl (y), (4.4) 

rj = ±1, and X = ex, Y = ey are long spatial scales of envelope functions A\ and A 2 . Substituting the above 
expansions into Eq. (|2.4[) . the equation at O(e) is automatically satisfied. At order 0(e 2 ), the equation for u\ is 

ui xx + uiyy - [F(x) + F(y)} m + Mo^i = -2 + ■ ( 4 - 5 ) 

Its homogeneous equation has two linearly independent solutions, pi(x)p 2 (y), and pi(y)p 2 (x) . In order for the inho- 
mogeneous equation (14. 5|) to admit a solution, the following Fredholm conditions 

rr(^ + s) piWKfa, ^ =o ' (4 - 6) 

must be satisfied. Here the integration length is 2L rather than L since the homogeneous solutions Pi(x)p 2 (y) and 
P2{x)pi(y) may have periods 2L along the x and y directions (see Eq. I|4.1j) ). Recalling the uq solution (|4.4|) . it is 
easy to check that the above Fredholm conditions arc indeed satisfied automatically, thus we can find a solution for 
Eq. (jHU) as 

9Ai , . . . dAi dA 2 dA 2 

u i = -QxM x m{y) + -QxMy)Pi{ x ) + -^-^wpitj/) + -^viiymw, (4.8) 

where v n (x) is a periodic solution of equation 

V n ,xx ~ F(x)v n + Uo^ n l> n = ~2p n ^ x , 71 = 1, 2. (4.9) 

At 0(e 3 ), the equation for u 2 is 

rn , . / 9 2 ui d 2 ui d 2 u d 2 u o \ ,. ..... 

w 2ra + w 2aa - [f (a) + F(y)] u 2 + Mo«2 = - I 2 g^r + 2 g^7 + fljp + ~g^2 + W + W M ° J ■ ( 4 - 10 ) 
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Substituting the formulas (|4.4[) and (|4.8[) for uq and u\ into this equation, we get 



{u 2x x + u 2yy - [F(x) + F(y)] u 2 + /X M 2 } 



|^ PM(z) +Pi(a:)]p2(y) + [2^(y) +»(»)] Pi (i) 
[2 ^ (x) +^W]pi(y) + ^ [2^(2/) +Pi(tf)]P2(a:) 



+p\(x) P l{y)\A 1 \ 2 A 1 + Pl (x)p 2 (y)p 2 1 (y)p 2 2 (x)(A 1 Al + 2^1 1 |A 2 | 2 ) 
+p3 ( X )pl ( y ) | A 2 1 2 A 2 + p? (x)pj (y )p 1 (y )p 2 (x) (A 2 A\ + 2 A 2 \ A 1 \ 2 ) 

+ v A 1 (X,Y) Pl (x) P2 (y) + r 1 A 2 (X,Y)p 2 (x)p 1 (y). (4.11) 

Here the overbar represents complex conjugation. Before applying the Fredholm conditions to this inhomogeneous 
equation, we notice the following identities: 



L 



2L 

Pl (x)p 2 (x)dx = 0, (4.12) 



2L p2L 

Pi{x)v 2 (x)dx — \ p 2 {x)vi{x)dx, (4-13) 
o Jo 



and 

r 2L r 2L 



[2v' n {x) +p n {x)]p n (x)dx = D n / pi(x)dx, n=l,2, (4.14) 
o Jo 



where 

_ld 2 u 
n= 2 



(4.15) 



Identity (|4.12p holds since p\{x) and p 2 (a;) are the eigenfunctions of the self-adjoint linear Schrodinger operator with 
different eigenvalues. Identity (|4. 13|) can be confirmed by taking the inner product between Eq. f|4. 0|) and functions 
p n (x). Identity (14.14[) can be verified by expanding the solution of Eq. (|3.3[) around the edge of the Bloch band 
to = wo,n (see Eq. (15) in [20]). Utilizing these identities and (14. ip . the Fredholm conditions for Eq. (|4.11|) finally 
lead to the following coupled nonlinear equations for the envelope functions A\ and A 2 : 

+ ° 2 W^ +^Ai+a [a\A 1 \ 2 A 1 + (3 (A 1 A 2 2 + 2yl 1 |A 2 | 2 ) + 7 (\A 2 \ 2 A 2 + A 2 A\ + 2A 2 \A 1 \ 2 )] = 0, (4.16) 
pa 



d 2 A 2 n d 2 A 2 
U2 f3x ' 

Here 



° 2 ~dX* + Dl W^ + VM + ° H^ 2 ^ 2 + 13 ( A * A * + 2 MM 2 ) + 7 (l^i| 2 ^i + MA 2 + 2A!|^ 2 | 2 )] - 0. (4.17) 



2L P 2L 



JO 



Pi ( x )P2(y) dxd v 



2L i-2L 



JO 



(4.18) 



Pi ( x )P2(v) <Jxdy 



p2L p2L 

/ / Pi( x )p 2 2(x)p 2 (y)pl(y) dxdy 

> = J ° \ 2 L f 2L > (4-19) 

/ p 2 ( x )pl(y) dxd v 

o Jo 
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and 

r 2L r 2L 



Pi (x)P2{x)p 2 (y)pi{y) dxdy 



2L r 2L 

2 



JO 



(4.20) 



Pi ( x )p 2 (y) dxdy 



Notice that a and are always positive, but 7 may be positive, negative, or zero. 

The coefficients in Eqs. (|4.16j) - (|4.17j) can be readily determined from Bloch solutions of the ID equation (|3.3p . In 
particular, when the ID Bloch waves pi{x) and P2{x) have been normalized to have unit amplitude (see Fig. 2), we 
find that at point C, 

Di = 0.434845, D 2 = 2.422196, a = 0.5821, = 0.1325, 7 = 0; (4.21) 

at point D, 

Di = -0.586799, D 2 = -13.264815, a = 0.5256, = 0.1811, 7 = 0; (4.22) 

and at point E, 

Di = 0.434845, D 2 = 15.793172, a = 0.4684, = 0.0781, 7 = -0.0261. (4.23) 

At band edges where a single Bloch mode exists (such as points A and B in Fig. [3fb)), the envelope equation for 
this single Bloch mode can be more easily derived. In this case, this single Bloch mode must be of the form pi (x)p\ (y), 
where pi(x) = p(x;u>o t i), and loq.i is a band edge in the ID problem (|3.3p . Unlike the previous case, interchanging x 
and y of this Bloch solution does not give different Bloch modes. The leading order solution 1*0(2;, y) in the expansion 
(|4.2p now is Ai(X,Y)pi(x)pi(y), and the envelope equation for A\[X, Y) can be found to be 

D 1 (^ + ^\+r ] A x +aa Q \Ai\*Ai=Q, (4.24) 

where Di is as given in Eq. (|4.15p . and 

r-2L r-2L 

\ \ Pi(x)pt(y) dxdy 

«o = J ~% J % • (4-25) 

Pi( x )Pi(y) dxdy 
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In particular, when the Bloch solution pi{x) is normalized to have unit amplitude, then at point A, 

Di = 0.434845, a = 0.462815; (4.26) 

and at point B, 

Di = -0.588073, a = 0.500922. (4.27) 
We note that for this single-Bloch-mode case, an envelope equation similar to (|4.24[) has been derived before in [37[. 



B. Locations of envelope solitons 



The envelope equations (|4.16p ~ (|4.17p and (|4.24p are translation-invariant. For instance, if [Ai(X,Y), A 2 (X,Y)] is 
a solution of Eqs. (|4~T6|) - (|4T7|) . so is [A X {X - X , Y - Y ),A 2 (X - X ,Y - Y )] where (X , Y ) are any constants. 
However, only when (Xo,Yq) take some special values can the perturbation series solution (|4.2p truly satisfy the 
original equation (|2.4p . The reason is that (Xq,Yq) must satisfy certain additional constraints. These constraints 
are exponentially small in e, thus they can not be captured in the power series expansions of (|4.2p . but need to be 
calculated using asymptotics beyond all orders techniques [HI or other equivalent methods [2(|. In ID problems, it 
has been shown that envelope solitons can only be located at two positions relative to the periodic lattice [2(| (or the 
underlying periodic wave train [38j]). In the present 2D problem, we will show below that envelope solitons can only 
be located at four positions relative to the 2D periodic lattice by a method similar to that used in [2p| . As we have 
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done in the previous subsection, detailed derivations will be presented for the {Ax, A 2 ) solutions of Eqs. I|4.16[) - (|4.17[) . 
Similar results hold for the A\ solution of Eq. (|4.24[) as well. 

First we derive two constraints for the envelope solutions. Multiplying Eq (|2.4| by u x or u y , adding its conjugate 
equation, and integrating from — 00 to +00, we get the following two constraints: 



-00 /"+oo 

|2 



00 «/ — 00 

+ OO Z' + OO 



— 00 J — 00 



F'{x)\u{x,y)\ 2 dxdy = 0, (4.28) 
F'(y)\u(x,y)\ 2 dxdy = Q. (4.29) 



Substituting the perturbation expansion (|4.2[) of the solution u{x,y) into the above equations, these constraints at 
the leading order become 



/ + 00 r + oc 

/ F'{x)\AxPx{x)p 2 {y) + A 2P2 (x) Pl (y)\ 2 dxdy = 0, (4.30) 
OO J — OO 

h{x ,y a ) = e 2 / F'{y)\AxPx{x)p 2 {y)+A 2 p 2 {x) P x{y)\ 2 dxdy = 0. (4.31) 



00 •/ — oc 



Here 

A k = A k (X -X a ,Y -Y Q ), k = l,2, (4.32) 

and (Xo,Y<j) — (exo,eyo) is the center position of the envelope solution (Ai,A 2 ). In this paper, we consider envelope 
functions (Ai,A 2 ) such that \Ai\ 2 , \A 2 \ 2 and A\A 2 + A\A 2 are symmetric in X and Y about the center position 
(X ,Y Q ), he. 

\A k (-X,Y)\ 2 = \A k (X,Y)\ 2 = \A k (X,-Y)\ 2 , k ~ 1,2, (4.33) 

and similar relations hold for AiA 2 + A\A 2 . All solitary- wave solutions of Eqs. (|4.16|) - (|4.17p that we know satisfy 
these assumptions (up to spatial translations). 

Now we examine the constraint (|4.30[) . This constraint can be rewritten as 

/ + OO <- + OC 
/ F'{x) [\Ax\ 2 p\{x)p 2 {y) + \A 2 \ 2 p 2 {x)p 2 {y) 
-OO J —OO 

+ {A X A 2 + AxA 2 )px{x)p 2 {x)px{y)p 2 {y)] dxdy = 0. (4.34) 

Since F'(x) is anti-symmetric, and p\ (x) , p\ (x) both symmetric, functions F' (x)p 2 (x)p 2 (y), F 1 (x)p\{x)p\{y) and 
F' {x)pi{x)p 2 {x)p\{y)p 2 {y) have the following Fourier series expansions: 



F'{ X )p\{x)p 2 {y) =J2J2 c«„ sin ^ cos (4.35) 

m— 1 n— 



F'{x)p 2 {x)p\{y) = £ £ ti] n sin cos (4.36) 

F'{x)px{x)p 2 {x)px{y)p 2 {y) = jr jr c% sin cos + £ £ cos sin ( 4 ' 37 ) 

m— 1 n=0 m— n=l 

("3") fS") 

Here dm,n = or Cm\ n — (for all m, n) if Pi(x)p2(^) is even or odd respectively. When the above Fourier expansions 
are substituted into Eq. (|4.34p , every Fourier mode in these expansions leads to an exponentially small term in (|4.34[) , 
and the exponential rate of decay of these terms is larger for higher values of m 4- n. Keeping only the leading-order 
term obtained from Fourier modes with m + n = 1, Eq. (|4.34[) becomes 



■OO Z' + OG 



OO •/ — OO 



h{x ,yo) = e' / { c\^\AxY+c^>\A 2 \' + c^'{AxA 2 +AxA 2 ) 



. 2nx 
sin — — 



+d { H{A 1 A 2 + A X A 2 ) sin ^ J dxdy. (4.38) 
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Recalling Eq. (|4.32[) and our symmetry assumptions on |^4i| 2 , \A 2 \ 2 and A\A 2 + A1A2 (see (|4.33p ). the above integral 
can be simplified to be 

h(x ,yo) = Wi,ism(2Trx /L)+Wi, 2 sm(2iry /L), (4.39) 

where 

/ + OO /*+0O - r\ 
I { c$ \A 1 (X,Y)\ 2 + eg \A 2 (X,Y)\ 2 + eg [A, (X, Y)A 2 (X,Y) + A 1 (X, Y)A 2 (X,Y)]\ cos -^dxdy , 

(4.40) 

and 

/ dfl [A 1 (X,Y)A 2 (X,Y)+A 1 (X,Y)A 2 (X,Y)] cos— ^dxdy. (4.41) 

-00 J —00 

Notice that both integrals VFi,2 and Wi,2 are exponentially small in e, thus the constraint (|4.30|) is exponentially 
small, hence it can not be captured by power-series perturbation expansions (|4.2p . Repeating similar calculations for 
the integral of I 2 (xo,yo) in Eq- l|4.3ip . we can get (to the leading order) 

h(x , Vo) = W 2 ,i sin(2wx /L) + W 2 , 2 sm(2ny /L), (4.42) 

where expressions for W 2 1 and W 22 are similar to those for W^i and Wi i2 above. Then in order for the two constraints 
(1430)) and P~3Tj) to hold, we must have 

sin(27ra;o/£) = sin(2vry /L) = 0. (4.43) 
Thus, the envelope solution (Ai,A 2 ) can only be centered at four locations: 

(x o ,y o ) = (0,0), (0,~), (|,0), (~ |). (4.44) 

Note that due to the x and y symmetry of the lattice, envelope solutions centered at the second and third locations 
as above are topologically the same (except an interchange of x and y axes). Other possible locations of (xq, yo), such 
as (i, L), are equivalent to one of the four locations above due to periodicity of the lattice potential. 



V. SOLUTIONS OF ENVELOPE EQUATIONS 

Envelope equations (14. 16j) - (14. 17[) and (|4.24p are one of the key results of this article. They have important conse- 
quences. First, they show that solitary waves are possible only when rjDi < and r\D 2 < 0. This simply means that 
\x must lie in the bandgap of the linear system (see Eq. (|4.3p V Second, these equations show that solitary waves exist 
only when the dispersion coefficients Di,D 2 and the nonlinearity coefficient a are of the same sign. For instance, 
at points A, C and E in Fig. [3] where D\ > and D 2 > 0, solitary waves exist only when a > 0, i.e., for focusing 
nonlinearity, but not for defocusing nonlinearity (cr < 0). The situation is opposite at points B and D. This result is 
consistent with our physical intuitions as well as the 2D experimental observations in [3l| . 

Below, we consider soliton solutions of envelope equations (|4.24|) and (|4. 16[1 - (|4. 1T|) . The scalar equation (14.241) is 
the familiar 2D NLS equation, and it admits the following types of solitary wave solutions: 

(1) A\ = f(R), which is real and radially symmetric; here (R, O) is the polar coordinates of (X, Y): 

(2) A\ — f(R)e l , where n is an integer; this is a vortex-ring solution with charge n. 

The coupled envelope equations (|4.16[) - (|4.17|) admit a wider array of solutions. A complete classification of their 
solutions is beyond the scope of the present article. Below we only list a few types of their solutions. If 7 = 0, these 
equations admit the following three simple solution reductions: 

(i) A\ 7^ 0, A 2 ~ 0; or A\ ~ 0, A 2 ^ 0. In this case, the solution is a single Bloch-wave envelope solution. Notice 
that A\ (or A 2 ) satisfies a 2D NLS equation with different dispersion coefficients D\ and D 2 along the X and Y 
directions, thus these solutions are simply those described above (charge- free or vortex solutions) except being 
stretched along the X or Y direction, hence becoming ellipse-shaped. 
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(ii) Ai, A2 £ R, where R is the set of real numbers. In this case, the solution is a real- valued vector envelope state. 
Note that if {A 1 ,A 2 ) £ R is a solution, so are (-At, A 2 ), {A x , -A 2 ), (-At, -A 2 ), and (iA 1 ,iA 2 ). All these 
solutions are equivalent to each other and lead to the equivalent solitary waves in the original system (|2.4p . 

(iii) Ai £ R, A 2 £ iR. In this case, the solution is a complex- valued vector envelope state. Note that the other 
solution of Ai £ iR, A 2 £ R is equivalent to this A\ £ R, A 2 £ iR solution. 

If 7 7^ 0, however, the reductions are quite different. For instance, the first and third reductions of case 7 = no 
longer hold. In this case, the reduction of A\, A 2 £ R is allowed. Under this reduction, there are two subcases, A\ > 0, 
A 2 > 0, and A\ > 0, A 2 < 0, which are not equivalent to each other. They lead to different solitary waves in Eq. 
(|2.4j) (see the end of this section) . 

It is note-worthy that at band edges where 7 7^ 0, the single Bloch-wave envelope reductions of A\ 7^ 0, A 2 = 
and Ai = 0, A 2 7^ are not possible. Physically, this is due to a resonance between the two Bloch modes, which 
prevents the existence of a single Bloch mode envelope solution. For instance, at point E of Fig. [3] where 7 7^ (see 
Eq. (|4.23p ). the two Bloch modes are both located at the same T point in the first Brillouin zone (see Sec. IIIII and 
Fig. StE)). These two Bloch modes are thus in resonance, which makes 7 7^ 0. At points where 7 = (such as points 
C and D in Fig. [31), the two Bloch solutions are located at different symmetry points of the Brillouin zone and are 
not in resonance, thus single Bloch-wave reductions of A\ 7^ 0, A 2 = and A\ — 0, A 2 7^ are possible there. 

To illustrate various envelope-soliton solutions admitted by Eqs. (|4.24|) and (|4.16[) - (|4.17p . we consider points 
A,B,C,D and E of Fig. [3] in detail below. At each of these five points, we determine solutions of the underlying 
envelope equations. For each envelope solution, we also display the corresponding leading-order analytical solution 
uo(x,y). Since envelope solitons can have four different locations (see the previous section), which will lead to four 
different solutions uq(x, y), for simplicity, we will only display the one where the envelope is centered at (xo, yo) = (0, 0) 
below. 

First we consider point A. At this point, a single Bloch mode exists and has been shown in Fig. IHA). This Bloch 
wave is located at T point of the first Brillouin zone, is 7r-periodic along both x and y directions, and is all positive. The 
envelope equation at this point is given by (|4.24[) with the coefficients D\ and ao given in Eq. (|4.26|) . Since D\ > 0, 
solitary waves will bifurcate into the semi- infinite bandgap (ry < 0) under focusing nonlinearity (a — I). If we take the 
solution of (|4.24[) to be A\ — f(R) where f(R) > is the ground-state envelope solution, the corresponding leading- 
order analytical solution uo(x,y) is a nodeless solitary wave. Such solutions have been reported before [1,0 
Other envelope solutions of (|4.24p such as A\ = f(R) where f(R) is a higher-mode solution (with nodes) or vortex-ring 
solutions Ai = /(i?)e me could lead to other types of solitary wave structures. 

Next, we consider point B. At this point, a single Bloch mode exists and has been shown in Fig. BJB). This Bloch 
mode is located at M point of the first Brillouin zone, is 27r-periodic along both x and y directions, and its adjacent 
peaks are out of phase. The envelope equation at this point is given by (|4.24p with the coefficients D\ and ao given 
in Eq. (|4.27p . Since D\ < 0, solitary waves will bifurcate into the first bandgap (77 > 0) under defocusing nonlinearity 
((j = — 1) . If we take the solution of (|4.24p to be the ground-state envelope solution, the corresponding analytical 
solution uo(x, y) is a gap soliton with nodes. Such solutions have been reported in [El. flBL [23l] . 

Envelope solutions at points C, D and E are richer and more interesting. Two of them, which are the so-called 
reduced-symmetry solitons [30j and gap vortex solitons [27) . have been reported before. But many other solutions 
at these points are novel and have not been considered yet. They include solutions which we call dipole- array 
solitons, dipole-cell solitons, vortex-cell solitons, etc. Envelope solutions and the corresponding leading-order analytical 
solutions uo(x,y) at these three points will be described in the next three subsections respectively. 

A. Envelope solutions at point C 

At point C, two Bloch modes exist. One of them is shown in Fig. |4jC), while the other is a 90°-rotation of Fig. 
0|C). These Bloch modes are located at the X and X' points of the first Brillouin zone. They are 7r-periodic along 
one spatial direction, and 27r-periodic along the other spatial direction. Envelope equations at this point are given by 
(|4~T61) - (|4~T7| . with the coefficients Di,D 2 ,a 7 (3 given in Eq. (|4.2ip . and 7 = 0. Since D\, D 2 > here, solitary waves 
will bifurcate into the first bandgap (rj < 0) under focusing nonlinearity (<r = 1). Since 7 = 0, we have three solution 
reductions (see above). Under these reductions, we consider the following subclasses of solutions. 

(i) A\ > 0, A 2 = 0. In this case, the envelope solution is an ellipse as shown in Fig. E)(a). This ellipse stretches 
along the Y direction. When the Bloch wave (see Fig. 0JC)) is modulated by this envelope function, the 
resulting leading-order analytical solution uo(x,y) is plotted in Fig. [D(b) (with e taken as 0.2). This solution 
contains only a single Bloch mode (since A 2 =0), thus we can call it a single- Bloch-mode soliton. This soliton 
is narrower along the x direction, and broader along the y direction, thus it is expected to be more mobile along 
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FIG. 5: (color online) (a) Envelope solution Ai(X, Y) at point C with A2 = 0; (b) the corresponding analytical single-Bloch- 
mode soliton uo(x,y) with e = 0.2. 



the y direction and less so along the x direction. This type of solution has been observed in [3(| for a different 
(saturable) nonlinearity (see Sec. I VII (A) for more details). In that paper, these solutions were called reduced- 
symmetry solitons. One of their potential applications in optical routing and switching has been described in 

(ii) A\ > 0, A 2 > 0. In this case, the envelope solutions are both real and positive, and they are shown in Fig. 
[6l[a, b). They are both ellipse-shaped but stretched along opposite directions. The corresponding leading-order 
analytical solution uo(x,y) is plotted in Fig. \@[c, d) (with e taken as 0.2). The central region of this solitary 
wave is a dipole array aligned along the two diagonal directions, thus we will call this solitary wave a dipole- 
array gap soliton. The outer region of this soliton is aligned along the horizontal (x) and vertical (y) directions 
though. Note that this dipole-array soliton is quite different from dipole solitons reported before (see [TiJ H(| for 
instance): previous dipole solitons reside in the semi-infinite bandgap, and their peaks are at lattice sites; but 
the present dipole-array soliton resides in the first bandgap, and its peaks are off lattice sites. This dipole-array 
gap soliton arises due to a superposition of two modulated Bloch modes (see Eq. (I4.4|l ). and has never been 
reported before. 

(hi) A\ > 0, A2 — iA2, A2 > 0. In this case, the envelope of one Bloch wave is real, while that of the other Bloch 
wave purely imaginary. In the literature, such two Bloch modes are called to have tt/2 phase delay [H], 
Envelope functions A\ and A2 look very similar to A\ and A2 of Fig. [^a, b). Indeed, it is easy to see that A\ 
and A2 satisfy the same equations as A\ and A2 of the previous reduction (ii) except that the /3 coefficient is 
slightly different. The leading-order analytical solution uo(x,y) for these envelope solutions is displayed in Fig. 
|()][e, f). This soliton looks quite different from the previous one in Fig. [^c, d). The most significant difference is 
that, when winding around each lattice center (i.e. points x = mit, y = rnr with m, n being integers), the phase 
of the present soliton increases or decreases by 2ir. In other words, the solution around each lattice center has a 
vortex structure. Thus we call this solution a vortex-array gap soliton. This vortex-array soliton is qualitatively 
the same as the gap vortex soliton reported in [27| for a different nonlinearity (see Sec. IVIf A)). However, it is 
quite different from the other vortex solitons residing inside the semi-infinite bandgap, as has been reported in 
[23, HI Hi HI before. 

For a vortex-type soliton, its angular momentum is an important quantity. This momentum is defined as 

//•oo />oo 
r x Im(fiVu) dr = Im / / u{xu y — yu x )dxdy . (5-1) 
J — 00 J — 00 

Here Im(uVit) is the linear momentum, and r is the position vector in the {x 1 y) plane. The spin of the vortex 
is defined as 

(5.2) 



In the absence of the periodic potential, the spin of a vortex soliton takes an integer value equal to the phase 
winding number (topological charge). In the present case with a periodic potential, the spin of a vortex-array 
gap soliton does not take integer values in general. Even though this soliton has a local vortex structure around 
each lattice center, the topological charges of adjacent lattice centers are opposite, thus the total spin of this 
vortex-array soliton is finite, not infinite. In fact, substituting the perturbation-series solution (|4.2p . (|4.4[) . and 
A\ e R, A 2 G iM. into the f2 and S formulas, we can easily show that both and S approach zero as e — > 0. 
Thus the spin of these vortex-array solitons is quite small. 
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FIG. 6: (color online) (a, b) Envelope solutions Ai(X, Y) > (left) and A2(X, Y) > (right) at point C; (c, d) the corresponding 
analytical dipole-array soliton Uo(x,y) with e = 0.2: the left is the amplitude and right the phase; (e, f) the analytical vortex- 
array soliton uo(x,y) at point C in the case (iii) with e = 0.2: the left is the amplitude and right the phase. 

B. Envelope solutions at point D 

Point D admits two Bloch waves, one shown in Fig. BJD), and the other one being a 90°-rotation of Fig. 0JD). 
These Bloch modes are located at the X and X' points of the first Brillouin zone. Envelope equations at this point 
are given by (|4.16|) - (|4.17p . with the coefficients Di,D 2 ,a,(3 given in Eq. (|4.22[) . and 7 = 0. Since D\,D 2 < now, 
solitary waves will bifurcate into the second bandgap (77 > 0) under defocusing nonlinearity {a = —1). Like point C, 
three solution reductions are admitted. Under these reductions, we consider the following subclasses of solutions. 

(i) Ai > 0, A 2 = 0. In this case, the envelope solution is an ellipse as shown in Fig. [7fa), which is thinner 
than Fig. 02a) at point C. When the Bloch wave (see Fig. 0jD)) is modulated by this envelope function, the 
resulting leading-order analytical solution uo(x,y) is plotted in Fig. EJb) (with e = 0.2). This solution is also a 
single-Bloch-mode soliton since A 2 — 0. Due to the thin envelope solution, this soliton is much broader in the 
y direction than in the x direction. Thus we can expect it to be much more mobile along the y direction than 
along the x direction. Note that this solution is under defocusing nonlinearity, while a similar-looking solution 
in Fig. [5(b) ( see a l so [30j |) was under focusing nonlinearity. 

(ii) A\ > 0, A 2 > 0. In this case, the envelope solutions are both real and positive, and are shown in Fig. [Ha, b). 
They are both ellipse-shaped stretching along opposite directions. The corresponding leading-order analytical 
solution uq(x, y) is plotted in Fig. [5fc, d) (with e = 0.2). This solution looks quite different from its counterpart 

- the dipole-array soliton of Fig. [5(c, d) at point C. Its main feature is that the solution inside each lattice 
site is a dipole cell. Its difference from the dipole-array soliton at point C is that, here the two humps of each 
dipole are completely confined inside each individual lattice, while the two humps of each dipole in Fig. [5Jc, d) 
spread to neighboring lattice sites. We will call this soliton in Fig. [5Jc, d) a dipole-cell gap soliton. Solutions of 
this kind have not been reported before. Notice that in the central region of the soliton, dipole cells are aligned 
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FIG. 7: (color online) (a) Envelope solution Ai(X,Y) at point D with A2 = 0; (b) the corresponding analytical single-Bloch- 
mode soliton uo(x,y) with e — 0.2. 




FIG. 8: (color online) (a, b) Envelope solutions Ai(X, Y) > (left) and A2(X, Y) > (left) at point D; (c, d) the corresponding 
analytical dipole-cell soliton uo(x,y) with e = 0.2: the left is the amplitude and right the phase; (e, f) the analytical vortex-cell 
soliton uo(x,y) at point D in the case (hi) with e = 0.2: the left is the amplitude and right the phase. 



along diagonal directions, but in the outer region, dipole cells are dominated and aligned along the horizontal 
and vertical directions. 

(iii) Ai > 0, A 2 = iA 2 , A 2 > 0. In this case, the two Bloch modes have ir/2 phase delay, and envelope functions 
(Ai,A 2 ) are similar to (Ai 7 A 2 ) of Fig. [8ja, b). The leading-order analytical solution uq(x, y) with these envelope 
solutions is displayed in Fig. [5Je, f). The most significant feature of this solution is that, in the central region, 
the solution inside each lattice site is a vortex cell (with charge 1 or —1). Thus we call it a vortex-cell gap soliton. 
Charges of adjacent cells are opposite, hence the total angular momentum and spin of this soliton are both finite. 
When e — > 0, the angular momentum and spin approach zero. Notice that each vortex cell here is completely 
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FIG. 9: (color online) (a, b) Envelope solutions Ai(X,Y) > (left) and Ai(X, Y) > (right) at point E; (c) the corresponding 
analytical spike-array soliton uo(x,y) with e = 0.2; (d) the analytical quadrupole-array soliton uo(x,y) at point E in the case 
(ii) with e = 0.2. 



isolated and confined inside each individual lattice, which is quite different from the vortex-array soliton in Fig. 
[6][e, f), where neighboring vortices are connected together and do not have clear-cut boundaries between them. 
In addition, this vortex-cell soliton is under defocusing nonlinearity, while the vortex-array soliton of Fig. [BJe, 
f) was under focusing nonlinearity (see also (27j). Furthermore, this vortex-cell soliton resides in the second 
bandgap, while the vortex-array soliton of Fig. [Qe, f ) resides in the first bandgap. Gap vortex solitons under 
defocusing nonlinearity as reported in [IH [23 | reside also inside the first bandgap; they are the counterparts 
of vortex solitons in the semi-infinite bandgap as reported in [2l|, [2~3 |. and are fundamentally different from 
the vortex-cell solitons bifurcated from point D here. Although these vortex-cell gap solitons have never been 
studied before, they do resemble the linear vortex-cell defect modes in photonic lattices as reported in (T^ |. 



C. Envelope solutions at point E 

Point E also admits two Bloch modes, one shown in Fig. 0jE), and the other being a 90°-rotation of Fig. 0|E). These 
Bloch modes are both located at T point of the first Brillouin zone, thus are 7r-periodic along both spatial directions. 
Envelope equations at this point are given by (|4. 16[) - (|4. . with the coefficients D\. D 2 , ct, (3, 7 given in Eq. (|4.23l) . 
Since D±, D2 > 0, solitary waves will bifurcate into the second bandgap (77 < 0) under focusing nonlinearity (a = 1). 
Since 7 ^ 0, the reduction of A±, A2 € K is allowed, and we consider the following two subclasses of solutions. 

(i) Ai > 0, A2 > 0. In this case, the envelope solutions are both real and positive, and they are shown in Fig. 
IHJa, b). They are both very thin ellipses stretched along opposite directions. The corresponding leading-order 
analytical solution uq(x, y) is plotted in Fig. HJc) (with e taken as 0.2). This solution has a pronounced 
cross-shape overall structure. In the central region of the (x, y) plane, the solution at each lattice center is an 
almost-circular positive spike, which recedes to negative backgrounds away from the lattice center. Thus we 
may call this solution a spike-array gap soliton. 

(ii) A\ > 0, A2 < 0. In this case, envelope functions Ai and | ^4.2 1 are similar to those of Fig. [i^a, b), thus not shown. 
The corresponding leading-order analytical solution uq(x, y) is plotted in Fig. [j^d) (with e = 0.2). This solution 
also has a pronounced cross shape. In its central region, the solution resembles an array of quadrupoles, thus 
we can call this solution a quadrupole-array gap soliton. 
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VI. FAMILIES OF SOLITARY WAVES BIFURCATED FROM BAND EDGES 



The above multi-scale perturbation analysis predicts various types of low-amplitude solitary waves when the prop- 
agation constant /i is near edges of Bloch bands. As the propagation constant moves away from band edges, these 
solitary waves will become more localized, and their amplitudes will increase. From an experimental point of view, 
more localized solitary waves are often easier to observe. For more localized solutions, the above perturbation analysis 
starts to break down, and numerical methods need to be employed. In this section, we compute whole families of soli- 
tary waves bifurcated from edges of Bloch bands. The numerical method we will use is the modified squared-operator 
iteration method developed in (35j . 

For illustration purpose, we determine families of solutions bifurcated from points C, D and E of Fig. [3fb). Leading- 
order analytical approximations on low-amplitude solutions of these families have been presented in Figs. [5]to[9] Note 
that solution families bifurcated from points A and B have been reported before in [l5|, [U [23| , thus they will not be 
computed in this article. 



At point C, three types of solitary waves have been predicted in the previous section under focusing nonlinearity: 
single-Bloch-mode solitons, dipole-array solitons, and vortex-array solitons, see Figs. EJb),[6jc, d) and[6](e, f) respec- 
tively. They all exist in the first bandgap, bifurcated from the edge point C of the second Bloch band. Numerically, 
we have obtained the entire families of these three types of solutions. Their power curves are displayed in Fig. 1101 
Here the power is defined as 



All the three power curves are non-monotonic. They have non-zero minimum values inside the bandgap, below which 
solitary waves of the respective family do not exist. This contrasts the ID case where solitary waves exist at all power 
levels [20j . Of these three power curves, the one for the single-Bloch-mode family is the lowest. Thus single-Bloch- 
mode solitons take the least power amount to excite. Solitary waves of all three families at /i = 7.189 and 6.189, which 
are 0.04 and 1.04 below the edge point C (where fj,o — 7.229), are plotted in Fig. [TT] Note that the shorter separation 
fi — fiQ = —0.04 is chosen so that fi — fig = — e 2 , where e = 0.2 is the value we have used when plotting all analytical 
solutions u (x,y) in Figs. [5] and [5] The solitons at /i = 7.189, shown in Fig. ITlTa. c, f), are weakly localized and 
have low amplitudes, and they are close to the band edge C. These numerical solutions are almost identical to the 
analytical ones shown in Figs. E3b),[|)Jc) and[5](e) (the phase fields of solitons in Fig. [TlTc. f) are also almost identical 
to the analytical ones in Fig. [BJd, f) and thus not shown). Thus the numerical results corroborate the analytical 
theory of the previous section. More significant solutions in Fig. [TTJ are the three solitons at the other propagation 
constant \x — 6.189, which is deep inside the first bandgap. These solutions, displayed in Fig. [TTTb). (d, e) and (g, 
h), are strongly localized. Indeed, the soliton of the single-Bloch-mode family in Fig. CQJb) almost becomes a single 
dipole aligned along the y axis at the lattice site of origin (x,y) = (0,0); the soliton of the dipole-array family in 
Fig- fTTTd. e) almost becomes a single dipole aligned along the y = —x direction at the lattice site of origin; and 
the soliton of the vortex-array family in Fig. CQJg, h) almost becomes a single vortex at the lattice site of origin. 
The phase structures of these strongly localized solitons, however, still resemble those of weakly localized ones (see 
Figs. [6jd, f) and [TlTe. h)). This strong localization of solitons in these solution families makes them amiable for 
on-axis excitation in photorefractive crystals. Indeed, the soliton of vortex-array family in Fig. [TTT g. h) looks almost 
identical to the gap vortex soliton observed in [27] ] with on-axis excitation. It should be noted that reduced-symmetry 
solitons observed in [3(| also bifurcate from single Bloch modes at point C. However, such solitons observed in [3(| at 
high powers look a little different from that in Fig. CQJb): their solitons are confined to three lattice sites, with the 
middle-site peak intensity higher than those at the two neighboring sites; while the soliton of Fig. CQJb) is confined to 
two lattice sites with equal peak intensities. The reason for this difference is the following. As we have explained in 
Sec. HVf B). envelopes of Bloch modes can be centered at four locations (see Eq. (|4.44[) h For the single-Bloch-mode 
soliton family shown in Figs. [TD1 and [TTTa. b), the envelope was centered at the origin (xo,yo) = (0,0), which leads 
to a dipole structure under strong localization. If the envelope is centered at (xo,yo) = (0,L/2) = (0,7r/2) instead, 
solitons under strong localizations would have a central peak shouldered by two equal but lower intensity peaks just 
like those observed in [30(. Thus, solitons observed in [30( belong to the single-Bloch-mode soliton family of point C 
with the envelope centered at (0,7r/2) rather than (0,0). From the experimental point of view, dipole-array solitons 
as shown in Fig. [TTT c.d.e). as well as the strongly localized single-Bloch-mode soliton of dipole- type in Fig. [TTJb), 
have never been observed before, and they still await experimental demonstration. 



A. Solution families bifurcated from point C 




(6.1) 
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FIG. 10: (color online) Power diagrams of the three families of solitary waves bifurcated from the edge point C (under focusing 
nonlinearity) . Bottom curve: single-Bloch-mode branch (marked by letter V); middle curve: dipole-array branch (marked by 
letter 'd'); top curve: vortex-array branch (marked by letter V). The marked thick dot points are 0.04 and 1.04 below the 
band edge C. 




FIG. 11: (color online) Weakly and strongly localized solutions on the three soliton branches bifurcated from the edge point 
C in Fig. [10] Propagation constants of these solutions are fi — 6.189 and 7.189, which are marked by thick dots in Fig. 1101 
Top row: single Bloch-mode solitons; u(x,y) is shown; (a) /i = 7.189; (b) fi — 6.189. Middle and bottom rows: dipole-array 
and vortex-array solitons respectively; (c, f) amplitudes (\u(x,y)\) at /j, = 7.189; (d, e) amplitude and phase of the dipole-array 
soliton at /i = 6.189; (g, h) amplitude and phase of the vortex-array soliton at fi = 6.189. 
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4 6 8 D 10 

FIG. 12: (color online) Power diagrams of the three families of solitary waves bifurcated from the edge point D (under defocusing 
nonlinearity) . Bottom curve: single-Bloch-mode branch (marked by letter 's'); middle curve: dipole-cell branch (marked by 
letter 'd'); top curve: vortex-cell branch (marked by letter V). The marked thick dot points are 0.04 and 0.54 above the band 
edge D. 



B. Solution families bifurcated from point D 

At point D, three types of solitary waves were predicted in Sec. [V] under defocusing nonlinearity: single-Bloch-mode 
solitons, dipole-cell solitons, and vortex-cell solitons, see Figs. E[b),|5[c, d) and[S](e, f) respectively. They all exist 
in the second bandgap, bifurcated from the edge point D of the second Bloch band. Numerically, we have obtained 
the entire families of these three types of solutions. Their power curves are displayed in Fig. 1121 These power curves 
are also non-monotonic and have non-zero minimum values inside the bandgap. Solitary waves of these three families 
at /i = 9.121 and 9.621, which are 0.04 and 0.54 above the edge point D (where /io = 9.081), are plotted in Fig. Q]|] 
The solitons at fj, = 9.121, shown in Fig. [l"37 a. c, f), are weakly localized and have low amplitudes. They are located 
close to the band edge D. These numerical solutions are almost indistinguishable from the analytical ones shown in 
Figs. EJb), (5(c) and[5] (e) where e = 0.2 (the phase fields of solitons in Fig. [TBT c, f) are also almost indistinguishable 
from the analytical ones in Fig. [5Jd, f) and thus not shown), thus numerical and analytical solutions near the band 
edge D are in agreement. Solutions at fi = 9.621 away from the edge point D are more localized. These solutions are 
displayed in Fig. [T^b). (d, e) and (g, h) respectively. The single-Bloch-mode soliton in Fig. [Tffib) is confined to only 
one lattice site along the x direction, but occupies quite a few lattice sites along the y direction. Thus this soliton will 
be highly mobile along the y direction and strongly trapped along the x direction. The dipole-cell soliton in Fig. [T3ld, 
e) is largely confined to the single lattice site at the origin as a dipole aligned along the y = —x direction, with long 
tails stretching along the horizontal and vertical (i.e. lattice) directions in a cross pattern. The vortex-cell soliton in 
Fig. [TBT g. h) is also largely confined to the single lattice site at the origin in the form of a vortex-ring with charge — 1, 
with long tails along the horizontal and vertical axes as well. None of these solitary wave structures in Fig. [13] was 
theoretically predicted or experimentally observed before. These structures exist under defocusing nonlinearity. Such 
nonlinearity can be obtained in photorefractive cr yst als with a negative bias charge or in Bose-Einstein condensates 
for certain types of atoms such as 87 Rb and 23 Na |14| . 



C. Solution families bifurcated from point E 

At point E, two types of solitary waves were predicted in Sec. [V] under focusing nonlinearity: spike-array solitons 
and quadrupole-array solitons, see Fig. E^c, d). They exist in the second bandgap, bifurcated from the edge point E 
of the third Bloch band. Numerically, we have obtained families of these two types of solutions. Their power curves 
are displayed in Fig. 1141 These power curves are also non-monotonic and have non-zero minimum values inside the 
bandgap. Solitary waves of these families at fj, = 9.77 and 9.33, which are 0.04 and 0.48 above the edge point E (where 
/io = 9.81), are plotted in Fig. [121 The solitons at fi = 9.77, shown in Fig. [TBTa. c), are weakly localized and have low 
amplitudes. They are almost the same as the analytical solutions in Fig. [SJc, d) where e = 0.2, which is expected. 
The solutions at \i = 9.33 away from the edge point E are more localized. These solutions are displayed in Fig. [T5T b. 
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FIG. 13: (color online) Weakly and strongly localized solutions on the three soliton branches bifurcated from the edge point 
D in Fig. [T2] Propagation constants of these solutions are fi = 9.121 and 9.621, which are marked by thick dots in Fig. 1121 
Top row: single Bloch-mode solitons; u(x,y) is shown; (a) fj, = 9.121; (b) n — 9.621. Middle and bottom rows: dipole-cell and 
vortex-cell solitons respectively; (c, f) amplitudes (|u(x, j/)|) at /i = 9.121; (d, e) amplitude and phase of the dipole-cell soliton 
at /i = 9.621; (g, h) amplitude and phase of the vortex-cell soliton at /i = 9.621. 



d). The localized solution of the spike-array soliton family in Fig. [TSJ^b) is now confined to only a few lattice sites on 
the x and y axes, with a dominant circular spike in its center. The localized solution of the quadrupole-array soliton 
family in Fig. I15f d) is confined also to only a few lattice sites on the x and y axes, with a quadrupole in its center. 
These types of solitary waves are novel and have not been reported before. 



VII. TIME-DEPENDENT ENVELOPE EQUATIONS 

In the previous sections, our focus was to obtain solitary waves in Eq. (|2.ip . thus our envelope equations (|4.1G|) - 
(|4.17[) and (|4.24[) were time-independent. If we look beyond solitary waves, we can further study how slowly modulated 
low-amplitude Bloch-wave packets evolve with time under nonlinear effects. In such a study, envelope equations of 
Bloch modes will be time-dependent. Such time-dependent envelope equations are essential not only for the tracking 
of Bloch-wave packets' movement, but also for the stability analysis of stationary envelope solutions obtained in Sec. 
fVl This is analogous to the role the nonlinear Schrodinger equation played in the study of water wave packets and 
light pulses in nonlinear optics [40l . l4~fl | . Derivation of time-dependent envelope equations is a simple extension of the 
analysis in Sec. HVi and will be performed briefly below. As before, the derivation will be carried out for the slightly 
more complicated case where the solution is near band edges with two different Bloch modes. Results for the simpler 
case of the solution being near band edges with a single Bloch mode will be given without elaboration. 

Consider the evolution of low-amplitude Bloch-wave packets in Eq. (|2.ip near a band edge fj, — fio with two different 
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FIG. 14: (color online) Power diagrams of the two families of solitary waves bifurcated from the edge point E (under focusing 
nonlinearity) . Upper curve: the spike-array branch (marked by letter V); lower curve: the quadrupole-array branch (marked 
by letter 'q'). The marked thick dot points are 0.04 and 0.48 below the band edge E. 
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FIG. 15: (color online) Weakly and strongly localized solutions on the two soliton branches bifurcated from the edge point E 
in Fig. 1141 Upper row: spike-array solitons; lower row: quadrupole-array solitons. Propagation constants of these solutions are 
/i = 9.77 (left column) and 9.33 (right column), which are marked by thick dots in Fig. 1141 



Bloch modes p\ (x)f>2 (y) and p2{x)p\{y) . The solution can be expanded into the following perturbation series: 

U{x, y, t) = e-^ ot (eU + e 2 ^ +e 3 U 2 + ...), (7.1) 

where 



U = A 1 (X, Y, 7>i (x)p 2 (y)+A 2 (X, Y, T)p 2 (x) Pl (y), (7.2) 

e <C 1, X = ex, Y = ey are slow spatial variables, and T = e 2 t is the slow time variable. When this expansion is 
substituted into Eq. (|2.1|) , we find at 0(e 2 ) that the equation for U\ is simply Eq. (|4. 5|) with (k = 0, 1) replaced 
by Uk, and with an additional term idUi/dt added onto its left hand side. The solution for Ux is still given by Eq. 
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(|4.8[) . At 0(e 3 ), we find that the equation for U 2 is Eq. (|4.11[) with Uk (k = 0, 1, 2) replaced by Iff., with an additional 
term idU 2 /dt added onto its left hand side, and with -qu replaced by idU /dT. Requiring the U 2 solution not to 
grow linearly with time t (suppression of secular terms), the following time-dependent equations will be derived for 
the envelope functions Ai(X, Y, T) and A 2 (X, Y, T): 



r) A r) 2 A r) 2 A 

+ Dl ~dX^ + D2 W^ + a ^ Al ^ Al + P ( AlA i + 2 MM 2 ) + ^ {\M 2 A 2 + MA\ + 2A 2 \A 1 \ 2 )] = 0, (7.3) 
8 A 8 2 A 8 2 A 

l ~d¥ + D2 ~dX^ + Dl ~dY^ + ° H^| 2 A 2 + (3 [A 2 A\ + 2AM 2 ) + 7 {\Ai\ 2 A 1 + A X A\ + 2^1 1 |A 2 | 2 )] = 0. (7.4) 

Here all coefficients D\, D 2 , a, (3 and 7 are the same as those given in Eqs. (|4.15[) and (|4.18|1 - (|4.20|) . 

If the solution of Eq. (12. 1|) is a low-amplitude Bloch-wave packet near a band edge fj, = /xq where a single Bloch 
mode pi(x)pi(y) exists, the perturbation expansion for the solution U is still Eq. (|7.ip . except that 

Uo = A 1 (X,Y,T) Pl (x) Pl (y) (7.5) 

now. In this case, the equation for the envelope function Ai(X, Y, T) reduces to 

dA 1 „ fd 2 A 1 d 2 A 1 \ . . l2 . 

1 ~W + Dl {~dX^ + -dY^J + aao ^ A ' = °» ( 7 - 6 ) 

where the coefficients D\ and ao are given in Eqs. (|4.15[) and (|4.25|) . 

Using the above time-dependent envelope equations, one can study the time evolution of Bloch-wave packets. In 
addition, one can analyze linear and nonlinear stabilities of stationary Bloch-wave packet solutions (i.e. solitary waves 
in Eqs. (|7.3|) - (|7.4[) and (|7.6|) ). Such studies lie outside the scope of the present article. 



VIII. SUMMARY AND DISCUSSION 



In this paper, we systematically studied various families of solitary waves which are bifurcated from Bloch-band 
edges in the two-dimensional NLS equation with a periodic potential. Near the band edges, we analytically derived 
envelope equations for low-amplitude Bloch-wave packets. Based on these envelope equations, many novel types of 
solitary waves inside bandgaps were predicted, including dipole-array solitons, dipole-cell solitons, vortex-cell solitons, 
etc. Away from the band edges, solitary waves were traced numerically, and strongly localized solitary waves of 
respective families were obtained. These solitary waves have distinctive intensity and phase patterns such as cross- 
shape intensity distributions and vortex-cell pha se profiles. Certain solitary waves previously reported in 2D photonic 
lattices such as reduced-symmetry solitons [30j | and gap vortex solitons [27| were found to be special cases of our 
general classes of solutions. 

In this paper, the relatively simple and obvious types of solutions in envelope equations were determined, and their 
induced solitary wave families computed. Weather other solutions of envelope equations would generate additional 
families of solitary waves is still an open question. For instance, the single-envelope equation (|4.24|) also admits vortex- 
ring solutions of the type f(R)e mB . Whether such solutions would lead to new families of vortex- type lattice solitons 
remains to be seen. For the coupled envelope equations (|4.16[) - (|4.17[) . it is not even clear what kind of additional 
solutions they admit beyond the ones we presented in this article. Thus a classification of solutions in the coupled 
envelope equations (|4.16|) - (|4. 17[) is highly desirable. Linear and nonlinear stabilities of solitary waves reported in this 
paper is another open question which merits careful investigation. From the experimental point of view, the challenge 
is to experimentally demonstrate the new types of solitary waves obtained in this article, such as dipole-array solitons 
and vortex-cell solitons. These open questions are beyond the scope of the present article, and will be left for future 
studies. 
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